Chapter 6 Beta diversity

load("data/data.Rdata")

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  dplyr::select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

6.1 location

6.1.1 Richness diversity

betadisper(beta_q0n$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     5 0.51463 0.102926 11.076    999  0.001 ***
Residuals 86 0.79917 0.009293                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                1.0000e-03 1.0000e-03 1.0000e-03 1.0000e-03    0.007
Brazil    8.1559e-05            1.0500e-01 7.0000e-02 9.1300e-01    0.129
CaboVerde 6.2922e-04 1.1254e-01            1.0000e-03 2.4000e-02    0.934
Spain     1.7363e-12 7.3505e-02 3.5940e-05            3.2000e-02    0.001
Denmark   1.3206e-07 9.0556e-01 2.9720e-02 2.3526e-02               0.056
Malaysia  9.6128e-03 1.3730e-01 9.2726e-01 4.6570e-04 5.9273e-02         

adonis2(beta_q0n$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	6.57607	0.2442751	5.559605	0.001
Residual	86	20.34468	0.7557249	NA	NA
Total	91	26.92075	1.0000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata$location, perm = 999)

                   pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 2.0520609 7.280171 0.19528265   0.001      0.015   .
2     Aruba vs CaboVerde  1 1.0026264 3.213494 0.09975614   0.001      0.015   .
3       Aruba vs Denmark  1 1.9513181 7.018855 0.19486614   0.001      0.015   .
4      Aruba vs Malaysia  1 1.1650328 3.678697 0.11257172   0.001      0.015   .
5         Aruba vs Spain  1 2.0433431 8.167535 0.21974917   0.001      0.015   .
6    Brazil vs CaboVerde  1 2.1090004 9.164611 0.24013374   0.001      0.015   .
7      Brazil vs Denmark  1 0.5703370 2.907951 0.09113563   0.001      0.015   .
8     Brazil vs Malaysia  1 0.9283278 3.953419 0.11996993   0.001      0.015   .
9        Brazil vs Spain  1 0.6343957 3.769479 0.11503016   0.001      0.015   .
10  CaboVerde vs Denmark  1 1.8099512 8.070070 0.22373314   0.001      0.015   .
11 CaboVerde vs Malaysia  1 1.2143821 4.593887 0.14094322   0.001      0.015   .
12    CaboVerde vs Spain  1 1.8141430 9.281722 0.24896172   0.001      0.015   .
13   Denmark vs Malaysia  1 1.0124763 4.418610 0.13629855   0.001      0.015   .
14      Denmark vs Spain  1 0.5305346 3.310771 0.10573905   0.001      0.015   .
15     Malaysia vs Spain  1 0.8450695 4.218747 0.13094076   0.001      0.015   .

#pdf("figures/beta_q0_loca.pdf",width=9, height=5)
beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.2 Neutral diversity

betadisper(beta_q1n$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     5 0.2728 0.054560 3.7079    999  0.004 **
Residuals 86 1.2654 0.014715                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                7.0000e-03 1.0000e-03 1.0000e-03 1.0000e-03    0.001
Brazil    9.2383e-03            7.8500e-01 1.9100e-01 6.5200e-01    0.852
CaboVerde 7.7521e-04 7.7033e-01            2.8900e-01 8.6500e-01    0.914
Spain     1.6646e-05 2.2327e-01 3.0424e-01            3.6500e-01    0.274
Denmark   3.3254e-04 6.6322e-01 8.7642e-01 3.7291e-01               0.760
Malaysia  1.1422e-03 8.5989e-01 8.9763e-01 2.4754e-01 7.7435e-01         

adonis2(beta_q1n$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	5.709562	0.221304	4.88821	0.001
Residual	86	20.090068	0.778696	NA	NA
Total	91	25.799630	1.000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata$location, perm = 999)

                   pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 1.3739907 4.789154 0.13766227   0.001      0.015   .
2     Aruba vs CaboVerde  1 1.0636589 3.806728 0.11603497   0.001      0.015   .
3       Aruba vs Denmark  1 1.6037782 5.813780 0.16699651   0.001      0.015   .
4      Aruba vs Malaysia  1 1.2584980 4.466640 0.13346545   0.001      0.015   .
5         Aruba vs Spain  1 1.4382789 5.558030 0.16083181   0.001      0.015   .
6    Brazil vs CaboVerde  1 1.5901673 7.027154 0.19505161   0.001      0.015   .
7      Brazil vs Denmark  1 0.6954159 3.122217 0.09719804   0.003      0.045   .
8     Brazil vs Malaysia  1 0.7313843 3.199010 0.09935119   0.001      0.015   .
9        Brazil vs Spain  1 0.4007935 1.948927 0.06297236   0.033      0.495    
10  CaboVerde vs Denmark  1 1.8198218 8.556140 0.23405480   0.001      0.015   .
11 CaboVerde vs Malaysia  1 1.3171418 6.019861 0.17695136   0.001      0.015   .
12    CaboVerde vs Spain  1 1.5415454 7.905417 0.22017337   0.001      0.015   .
13   Denmark vs Malaysia  1 1.0691578 4.970172 0.15074753   0.001      0.015   .
14      Denmark vs Spain  1 0.5766622 3.014217 0.09718823   0.001      0.015   .
15     Malaysia vs Spain  1 0.6446972 3.265586 0.10444666   0.001      0.015   .

#pdf("figures/beta_q1n_loca.pdf",width=9, height=5)
beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.3 Phylogenetic diversity

betadisper(beta_q1p$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq Mean Sq      F N.Perm Pr(>F)  
Groups     5 0.2152 0.04304 2.1047    999  0.069 .
Residuals 86 1.7587 0.02045                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                0.06300000 0.06400000 0.00200000 0.10100000    0.002
Brazil    0.06236518            0.85700000 0.32700000 0.65800000    0.709
CaboVerde 0.05279004 0.86912737            0.18700000 0.73000000    0.537
Spain     0.00060645 0.32391452 0.20693940            0.10900000    0.500
Denmark   0.09811567 0.63811338 0.73976492 0.10108967               0.331
Malaysia  0.00519311 0.67845753 0.52055144 0.50696000 0.31212601         

adonis2(beta_q1p$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	2.690791	0.2683828	6.309563	0.001
Residual	86	7.335152	0.7316172	NA	NA
Total	91	10.025944	1.0000000	NA	NA

pairwise.adonis(beta_q1p$S, sample_metadata$location, perm = 999)

                   pairs Df  SumsOfSqs    F.Model         R2 p.value p.adjusted sig
1        Aruba vs Brazil  1 0.64714035  6.1133065 0.16928127   0.001      0.015   .
2     Aruba vs CaboVerde  1 0.56229778  5.3504131 0.15575979   0.001      0.015   .
3       Aruba vs Denmark  1 1.03758657  9.5591181 0.24790811   0.001      0.015   .
4      Aruba vs Malaysia  1 0.48047063  5.0352414 0.14794199   0.002      0.030   .
5         Aruba vs Spain  1 0.60047832  6.7919744 0.18976250   0.001      0.015   .
6    Brazil vs CaboVerde  1 0.70961443  8.1023431 0.21837821   0.001      0.015   .
7      Brazil vs Denmark  1 0.32742702  3.5968626 0.11034383   0.012      0.180    
8     Brazil vs Malaysia  1 0.19570506  2.5119802 0.07971509   0.027      0.405    
9        Brazil vs Spain  1 0.06456161  0.9106384 0.03044530   0.453      1.000    
10  CaboVerde vs Denmark  1 1.07120011 11.9405505 0.29895808   0.001      0.015   .
11 CaboVerde vs Malaysia  1 0.50032095  6.5728056 0.19011490   0.001      0.015   .
12    CaboVerde vs Spain  1 0.82104026 11.9236999 0.29866220   0.001      0.015   .
13   Denmark vs Malaysia  1 0.54258755  6.8084757 0.19559821   0.001      0.015   .
14      Denmark vs Spain  1 0.24626102  3.3999417 0.10827860   0.024      0.360    
15     Malaysia vs Spain  1 0.26932060  4.5771932 0.14050299   0.001      0.015   .

#pdf("figures/beta_q1p_loca.pdf",width=9, height=5)
beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.1.4 Functional diversity

betadisper(beta_q1f$S, sample_metadata$location) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     5 1.5016 0.300325 9.1732    999  0.001 ***
Residuals 86 2.8156 0.032739                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Aruba     Brazil  CaboVerde      Spain    Denmark Malaysia
Aruba                1.0000e-03 5.9400e-01 1.0000e-03 1.0000e-03    0.497
Brazil    8.6955e-04            2.0000e-03 6.8900e-01 2.2300e-01    0.002
CaboVerde 5.8936e-01 9.7727e-04            3.0000e-03 1.0000e-03    0.892
Spain     1.5125e-03 6.8358e-01 1.5804e-03            3.9000e-02    0.002
Denmark   1.7539e-04 2.3957e-01 9.2505e-05 3.3658e-02               0.001
Malaysia  4.9615e-01 1.0071e-03 8.8292e-01 1.5959e-03 7.4748e-05         

adonis2(beta_q1f$S ~ location, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
location	5	8.815599	0.5844961	24.19552	0.001
Residual	86	6.266793	0.4155039	NA	NA
Total	91	15.082391	1.0000000	NA	NA

pairwise.adonis(beta_q1f$S, sample_metadata$location, perm = 999)

                   pairs Df   SumsOfSqs    F.Model           R2 p.value p.adjusted sig
1        Aruba vs Brazil  1  3.77109355 40.6261637  0.575228238   0.001      0.015   .
2     Aruba vs CaboVerde  1  0.81599223  5.5801053  0.161367505   0.011      0.165    
3       Aruba vs Denmark  1  4.54372380 50.8948430  0.637022880   0.001      0.015   .
4      Aruba vs Malaysia  1  0.87672195  6.0397240  0.172367910   0.021      0.315    
5         Aruba vs Spain  1  3.49955224 36.9984243  0.560595570   0.001      0.015   .
6    Brazil vs CaboVerde  1  1.79573937 27.7696018  0.489163230   0.001      0.015   .
7      Brazil vs Denmark  1  0.27026868 35.0545409  0.547260825   0.001      0.015   .
8     Brazil vs Malaysia  1  1.69391038 26.6369446  0.478763612   0.001      0.015   .
9        Brazil vs Spain  1 -0.01258568 -0.9666596 -0.034482498   0.953      1.000    
10  CaboVerde vs Denmark  1  2.85656348 47.6207158  0.629731090   0.001      0.015   .
11 CaboVerde vs Malaysia  1 -0.02164970 -0.1836835 -0.006603446   0.986      1.000    
12    CaboVerde vs Spain  1  1.56212224 23.8545897  0.460028512   0.001      0.015   .
13   Denmark vs Malaysia  1  2.82610262 48.0023359  0.631590271   0.001      0.015   .
14      Denmark vs Spain  1  0.38029019 58.5480612  0.676480332   0.001      0.015   .
15     Malaysia vs Spain  1  1.47580731 22.9256154  0.450178466   0.001      0.015   .

#pdf("figures/beta_q1f_loca.pdf",width=9, height=5)
beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = "sample") %>%
  group_by(location) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = location, fill = location)) +
    geom_point(size = 4) +
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    scale_color_manual(values = location_colors)+
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

#dev.off()

6.2 Behaviour

6.2.1 Richness diversity

betadisper(beta_q0n$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00411 0.0041147 0.2442    999  0.635
Residuals 90 1.51634 0.0168483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.636
Feral     0.62238      

adonis2(beta_q0n$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.4407951	0.0163738	1.498173	0.176
Residual	90	26.4799587	0.9836262	NA	NA
Total	91	26.9207538	1.0000000	NA	NA

pairwise.adonis(beta_q0n$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.4407951 1.498173 0.0163738    0.09       0.09    

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

6.2.2 Neutral diversity

betadisper(beta_q1n$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0080 0.0080048 0.5581    999   0.47
Residuals 90 1.2909 0.0143435                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.474
Feral     0.45698      

adonis2(beta_q1n$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.386668	0.01498735	1.369385	0.315
Residual	90	25.412962	0.98501265	NA	NA
Total	91	25.799630	1.00000000	NA	NA

pairwise.adonis(beta_q1n$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1  0.386668 1.369385 0.01498735   0.165      0.165    

pdf("figures/beta_q1n_behaviour.pdf",width=9, height=5)
beta_q1n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 
dev.off()

quartz_off_screen 
                2 

6.2.3 Phylogenetic diversity

betadisper(beta_q1p$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.01229 0.012295 0.5442    999  0.471
Residuals 90 2.03337 0.022593                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic          0.474
Feral     0.46262      

adonis2(beta_q1p$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.2191073	0.02185403	2.010807	0.167
Residual	90	9.8068363	0.97814597	NA	NA
Total	91	10.0259435	1.00000000	NA	NA

pairwise.adonis(beta_q1p$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.2191073 2.010807 0.02185403   0.087      0.087    

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) 

6.2.4 Functional diversity

betadisper(beta_q1f$S, sample_metadata$origin) %>% permutest(., pairwise=TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.2486 0.248604 3.3502    999  0.072 .
Residuals 90 6.6786 0.074206                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Domestic Feral
Domestic           0.07
Feral    0.070508      

adonis2(beta_q1f$S ~ origin, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))), 
        permutations = 999,
        strata = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))) %>% pull(location)) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
origin	1	0.3816603	0.02530502	2.336579	0.173
Residual	90	14.7007312	0.97469498	NA	NA
Total	91	15.0823914	1.00000000	NA	NA

pairwise.adonis(beta_q1f$S, sample_metadata$origin, perm = 999)

              pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1 Domestic vs Feral  1 0.3816603 2.336579 0.02530502   0.113      0.113    

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(origin) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = origin, fill = origin)) +
    geom_point(size = 4) +
    scale_color_manual(values = origin_colors) +
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    )